##############################################################################################################################
## Fatality Thresholds, Causal Heterogeneity, and Civil War Research: Reconsidering the Link Between Narcotics and Conflict ##
##############################################################################################################################
################################################ Noel Anderson & Alec Worsnop ################################################
##############################################################################################################################
###################################################### Code for Figures ######################################################
##############################################################################################################################

## Install required packges
install.packages(c("foreign","cem","mvtnorm","MASS","lmtest","MatchIt",
                 "optmatch","sandwich","lme4","car","fields","stargazer","survival"))

# Note: Newer versions of the Zelig package no long support cox proportional hazard models. We have included
# an older version of Zelig (V. 3.5.4) in our replication materials on Dataverse. To install the older version 
# of Zelig, download the .tar.gz file and run the following:
install.packages("<<insert file location>>/Zelig_3.5.4.tar.gz", repos = NULL, type = "source")

## Load libraries
library(foreign)
library(cem)
library(Zelig)
library(mvtnorm)
library(MASS)
library(lmtest)
library(MatchIt)
library(optmatch)
library(sandwich)
library(lme4)
library(car)
library(fields)
library(stargazer)
library(survival)

## Load journal survey dataset
survey.data <- read.csv("<<insert file location>>/journalsurvey.csv")

## Load analysis dataset
dataset <- read.dta("<<insert file location>>/dataset.dta")

## Create necessary functions

# (1) Function to get robust variance covariance matrix from clustered SEs
clv <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  return(dfc*sandwich(fm, meat=crossprod(uj)/N))
}

# (2) Function to complete regressions with clustered standard errors and simulate percent change in dependent variable. 
# This function is composed of three elements: the independent variables, the dataset, and indicator for the treatment variable.  
# Note that the DV should be the first column in the data given to the function.

qsim <- function(X,data,treat){
  # Regress data and given covariates
  M2b <- lm(data[,1]~as.matrix(X), data=data) 
  # Generate a robust covariance matrix, clustered by country
  rvarcov <- clv(data,M2b,data$ccode)
  # Create design matrix
  Xdes <- as.matrix(cbind(1,X))
  # Create covariates values to simulate percent change
  covs <- as.matrix(apply(Xdes,2,mean))
  # Create "treatment"
  covs[treat,]<-1
  # Create "control"
  covs2 <- covs
  covs2[treat,] <- 0
  # Create empty object for the treated dataset
  simdata <- c()
  # Create empty object for the control dataset
  simdata2 <- c()
  for(i in 1:1000){
    # Simulate betas from the original model
    simbetas <- rmvnorm(1,mean=M2b$coefficients,sigma=rvarcov)
    simsigma <- sqrt(M2b$df.residual* summary(M2b)$sigma^2 / rchisq(1, M2b$df.residual))
    # Generate the simulated percent change values
    simdata[i] <- (mean(rnorm(n=1000,(simbetas%*%covs),(simsigma))))
    simdata2[i] <- (mean(rnorm(n=1000,(simbetas%*%covs2),(simsigma))))
  }
  # Create vector of percent change 
  TEs <- (exp(simdata)-exp(simdata2))/exp(simdata2)
  # Take average of the 1,000 percent change values
  ate <- (mean(TEs))
  # Confidence intervals
  low <- (quantile(TEs,.025))              
  high <- (quantile(TEs,.975))  
  low90 <- (quantile(TEs,.05))              
  high90 <- (quantile(TEs,.95)) 
  # Print the percent change and confidence intervals
  cbind(ate,low,high,low90,high90)
}

###########################################
##### FIGURE 1 A & B (Journal Survey) #####
###########################################

# Make indicator for whether an article used multiple thresholds numeric
survey.data$change=as.numeric(survey.data$change)
                    
# Generate yearly totals for high, low, and multiple thresholds (low, hlb=1; high, hlb=2; both hlb=3)
high=c()
low=c()
both=c()
change2=c()
for (i in 2000:2013){
  high=c(high,sum(survey.data$year==i & survey.data$hlb==2))
  low=c(low,sum(survey.data$year==i & survey.data$hlb==1))
  both=c(both,sum(survey.data$year==i & survey.data$hlb==3))
  change2=c(change2,sum(survey.data$year==i & survey.data$change==1,na.rm=T))
  }
                    
# Create totals for high, low and multiple
totallow=sum(low)
totalhigh=sum(high)
totalboth=sum(both)
totalchange=sum(change2)
total=totallow+totalhigh+totalboth
                    
# Create joint high/low category as well as only those that changed
highlow=c(totallow/total,totalhigh/total,0)
both=c(0,0,totalboth/total)
bardata=cbind(highlow,both)
                    
# Create bar plot
par(mar=c(5,5,4,2),mgp=c(2.2,1,0))
barplot(bardata,ylim=c(0,.8),axes=F,ylab="Percentage of Published Articles Employing\nLarge-N Civil War Datasets",
  names.arg=c("Only Low or High\nThreshold","High and Low\nThreshold"),
  col=c("gray46","gray56","gray80"),border=NA,space=.1)
  axis(2,seq(0,.8,.1),label=paste(seq(0,80,10),"%",sep=""))
  legend("topright",inset = c(.06, .1),  fill = c("gray46","gray56","gray80"),bty="n",cex=.8,border=NA,lty=0,
  legend = c("Only Low Threshold","Only High Threshold","High and Low Threshold"))
  text(.61,0.4938272/2, labels="n=77")
  text(.61, 0.4938272 + (0.2901235/2), labels="n=46")
  text(1.71, 0.2160494/2, labels="n=35")
                    
# Create pie chart
par(mar=c(5,5,4,2),mgp=c(2.2,1,0),xpd=T)
pie(c(totalchange/totalboth,1-totalchange/totalboth),angle=c(0, 160),
  labels=NA,clockwise=F,col=c("white","black"))
  text(0,.35,"51%",cex=1.5)
  text(0,-.35,"49%",col="white",cex=1.5)
  legend("bottom",c("Main results unchanged","Main results change"), cex=0.8, angle=c(0, 160),density=c(0,500),col=c("white","black"),lty=0,bty="n")
                   
##########################################
##### FIGURE 2 (Total Battle Deaths) #####
##########################################

# Select variables for the analysis and pare the dataset down
keep.vars <- c("BTTLDtotLN","ccode","DRUGP","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")
all.data <- dataset[keep.vars]
data <- na.omit(all.data)

# Select treatment
treat="DRUGP"

# Set thresholds from 25 to 1,000 in increments of 25
btldt=seq(25,1000,25)

# Create empty objects for the estimates of the first difference, the 90 and 95 Confidence intervals, the Number of obs, and the Number of treated obs
fd.est=c()
low.est=c()
high.est=c()
low90.est=c()
high90.est=c()
n.est=c()
t.est=c()

for (i in btldt){
  # Create dataset for the new battle death threshold (i.e. keep data with only the threshold or higher)
  datatest=data[exp(data$BTTLDtotLN)>i,]
  # Select independent variables
  X=datatest[,c("ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")]
  # Get percent change and confidence intervals going from control to treated by using the qsim function created above. 
  set.seed(4321)
  tester=qsim(X,datatest,treat)
  fd.est=c(fd.est,tester[1])
  low.est=c(low.est,tester[2])
  high.est=c(high.est,tester[3])
  low90.est=c(low90.est,tester[4])
  high90.est=c(high90.est,tester[5])
  n.est=c(n.est,dim(datatest)[1])
  t.est=c(t.est,sum(datatest$treat))
  print(i)
}

# Create the figure by plotting the mean percent changes and confidence intervals
par(xpd=FALSE)
plot(btldt,fd.est,ylim=c(min(low.est),max(high.est)),ylab="Percent Change in Battle Deaths \n(No Drug Cultivation in Conflict Zone to Drug Cultivation in Conflict Zone)"
     ,pch=20,xlab="Battle Death Threshold",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)
lines(btldt,low.est,lty=1,col="gray39",lwd=2)
lines(btldt,high.est,lty=1,col="gray39",lwd=2)
lines(btldt,low90.est,lty=1,col="gray",lwd=2)
lines(btldt,high90.est,lty=1,col="gray",lwd=2)
axis(1, at=seq(0,1000,50), labels=seq(0,1000,50),cex.axis=.7)
axis(2, at=seq(round(min(low.est),1),round(max(high.est),1),.2),cex.axis=.8,labels=seq(100*round(min(low.est),1),100*round(max(high.est),1),100*.2))
abline(h=0,lty=1)
legend(x=7,y=1.45,legend=c("95% Confidence Interval","90% Confidence Interval"),lty=c(1,1),lwd=c(1.8,1.8), col=c("gray39","gray"),cex=.5)

########################################################################################################
##### FIGURE 3: (Dropping four observations with drug cultivation and fewer than 50 battle deaths) #####
########################################################################################################

# Select variables for the analysis and pare the dataset down
keep.vars <- c("BTTLDtotLN","ccode","DRUGP","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")
all.data <- dataset[keep.vars]
data <- na.omit(all.data)

## Analysis with full data set
X=data[,c("ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")]
# Get percent change and confidence intervals going from control to treated by using the qsim function created above
set.seed(4321)
tester=qsim(X,data,treat)
fd.est=tester[1]
low.est=tester[2]
high.est=tester[3]
low90.est=tester[4]
high90.est=tester[5]

## Analysis with four observations with drug cultivation and fewer than 50 battle deaths dropped
# Drop observations with drug cultivation and fewer than 50 battle deaths
data4=data[!(exp(data$BTTLDtotLN)<50&data$DRUGP==1),]

X=data4[,c("ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")]
# Get percent change and confidence intervals going from control to treated by using the qsim function created above
set.seed(4321)
tester=qsim(X,data4,treat)
fd4.est=tester[1]
low4.est=tester[2]
high4.est=tester[3]
low490.est=tester[4]
high490.est=tester[5]

# Plot Figure 3
par(xpd=FALSE)
par(mar=c(5,5,2,2)) 
par(mgp=c(2.2,1,0))
plot(c(-.1,.1),c(fd.est,fd4.est),ylim=c(-.81,.41),ylab="Percent Change in Battle Deaths \n(No Drug Cultivation in Conflict Zone to Drug Cultivation in Conflict Zone)",
  pch=20,cex.lab=.8,cex.axis=.7,yaxt = "n",xaxt = "n",xlim=c(-.2,.2),type='n',xlab="")
arrows(-.1,low.est,-.1,high.est,lty=1,col="gray39",lwd=2.4, angle=90, code=3,length = 0.05)
arrows(.1,low4.est,.1,high4.est,lty=1,col="gray39",lwd=2.4, angle=90, code=3,length = 0.05)
arrows(-.1,low90.est,-.1,high90.est,lty=1,col="gray",lwd=2.4, angle=90, code=3,length = 0.05)
arrows(.1,low490.est,.1,high490.est,lty=1,col="gray",lwd=2.4, angle=90, code=3,length = 0.05)
points(c(-.1,.1),c(fd.est,fd4.est),pch=20,cex=1.2)
abline(h=0,lty=1)
legend(x=-.19,y=.25,legend=c("95% Confidence Interval    ","90% Confidence Interval"),lty=c(1,1),lwd=c(2,2), col=c("gray39","gray"),cex=.66)
axis(2, at=seq(round(min(low.est),1),round(max(high4.est),1),.1),labels=seq(100*round(min(low.est),1),100*round(max(high4.est),1),100*.1))
abline(h=0,lty=1)
mtext("Excluding the Conflicts with Drug\nCultivation and Fewer than 50 Battle Deaths\n(N=240)",line=1.804,side=1,at=.1,cex=.8)
mtext("Including All Observations\n(N=244)",line=1,side=1,at=-.1,cex=.8)

##############################################################
##### Figure 4 (Random Drop Test for Total Battle Deaths #####
##############################################################

# Here we simulate the effect of simply randomly deleting as many observations as are deleted at each 
# threshold above by sampling 1000 datasets for each threshold

# Select variables for the analysis and pare the dataset down
keep.vars <- c("BTTLDtotLN","ccode","DRUGP","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")
all.data <- dataset[keep.vars]
data <- na.omit(all.data)

# Get unique set of number of deleted observations at each threshold
index=c()
for (i in 2:length(n.est)-1){
  index[i]=(dim(data)[1]-n.est[i])
}
deld.est=unique(index)

# Create empty objects for the percent changes and confidence intervals
fdd.rans=c()
lowd.rans=c()
highd.rans=c()
lowd.rans90=c()
highd.rans90=c()
n.rans=c()
t.rans=c()

# Start from the 2nd observation as the first is zero
for (i in 2:length(deld.est)){
  # Create empty objects for each of the 1,000 datasets to be produced at each drop value
  fd.ransin=c()
  low.ransin=c()
  high.ransin=c()
  low.ransin90=c()
  high.ransin90=c()
  n.ransin=c()
  t.ransin=c()
  # Now create 1000 samples of deleted observations
  for (j in 1:1000){ 
    # Randomly select the number of data points to be deleted
    set.seed(j+(i*1000))
    rand=sample(dim(data)[1],deld.est[i], replace=F) 
    # Remove randomly selected observations from the dataset
    datatest=data[-rand,]
    X=datatest[,c("ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")]
    # Record percent change and confidence intervals for each dataset
    tester=qsim(X,datatest,treat)
    fd.ransin=c(fd.ransin,tester[1])
    low.ransin=c(low.ransin,tester[2])
    high.ransin=c(high.ransin,tester[3])
    low.ransin90=c(low.ransin90,tester[4])
    high.ransin90=c(high.ransin90,tester[5])
    n.ransin=c(n.ransin,dim(datatest)[1])
    t.ransin=c(t.ransin,sum(datatest$DRUGP))
  }
  fdd.rans=c(fdd.rans,mean(fd.ransin))
  lowd.rans=c(lowd.rans,mean(low.ransin))
  highd.rans=c(highd.rans,mean(high.ransin))
  lowd.rans90=c(lowd.rans90,mean(low.ransin90))
  highd.rans90=c(highd.rans90,mean(high.ransin90))
  n.rans=c(n.rans,mean(n.ransin))
  t.rans=c(t.rans,round(mean(t.ransin),0))
  print(i)
}

# Combine the battle death thresholds with the number of observations and treated units
summary=cbind(btldt,n.est,t.est)
summary=rbind(1,summary)
x=matrix(,dim(summary)[1],3)
for (i in 2:dim(summary)[1]){
  x[i,1:3]=ifelse(summary[i,2]==summary[i-1,2],NA,i)
}

# Get unique number of deleted observations
temporary=na.omit(x[,1])
uniquend=(summary[temporary,])

# Create the plot by graphing the mean percentage change across each of the 1000 datasets at each drop threshold
par(mar=c(5,5,4,2),mgp=c(2.5,1,0))
par(mgp=c(2.2,1,0))
plot(1:28,fdd.rans,ylim=c(-.84,1.42),
     ylab="Percent Change in Total Battle Deaths\n(No Drug Cultivation in Conflict Zone to Drug Cultivation in Conflict Zone)"
     ,pch=20,
     xlab="Number of Observations Randomly Deleted",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)
lines(1:28,lowd.rans,lty=1,col="gray39",lwd=2)
lines(1:28,highd.rans,lty=1,col="gray39",lwd=2)
lines(1:28,lowd.rans90,,lty=1,col="gray",lwd=2)
lines(1:28,highd.rans90,lty=1,col="gray",lwd=2)
axis(1, at=seq(1,28,1), labels=F,cex.axis=.7,tick=T)
axis(1, at=seq(1,28,2), labels=deld.est[ seq(2,28,2)],cex.axis=.7,tick=F)
axis(2, at=seq(-.8,1.4,.2),cex.axis=.7,labels=seq(-80,140,100*.2))
legend(x=.5,y=1.45,legend=c("95% Confidence Interval","90% Confidence Interval"), lty=c(1,1),lwd=c(1.8,1.8), col=c("gray39","gray"),cex=.66)
axis(3, at=1:28,label=F,cex.axis=.7)
axis(3, at=seq(1,28,2),label=uniquend[seq(2,28,2),1],cex.axis=.7,tick=F)
mtext("Corresponding Battle Death Threshold", side=3, cex=.8,line=2)
abline(h=0,lty=1)

########################################
##### FIGURE 5 (Conflict Duration) #####
########################################

# Create resources dummy for conflicts with hydrocarbons in and outside the conflict zone and with gem production
resources=c()
for (i in 1:dim(dataset)[1]){
  resourcevar=ifelse(dataset[i,"OUTZhydrocarbonP"]==1|dataset[i,"hydrocarbonP"]==1|dataset[i,"ALLGEMP"]==1,1,0)
  resources=c(resources,resourcevar)
}
dataset$resources=resources

# Select variables for the analysis and pare the dataset down
keep.vars <- c("BTTLDtotLN","BTTLDtot","resources","ongoingCONF","Duration","DRUGP","mountain","confbord","Weakrebel","Equalstrenght","coldwar",
               "internatz","popLN","ethfrac","DEM","ccode","terr","gdplLN")  
all.data <- dataset[keep.vars]
data <- na.omit(all.data)

# Create censored observations
data$ongoingCONF=1-data$ongoingCONF

# # Set thresholds from 25 to 1,000 in increments of 25
btldt=seq(25,1000,25)

# Create empty objects for percent change, confidence intervals, number of observations, and number of treated observations
fd.est=c()
low.est=c()
high.est=c()
low90.est=c()
high90.est=c()
t.est=c()
nd.est=c()

# Conduct analysis
for (i in btldt){
  # Create dataset for the new battle death threshold (i.e. keep data with only the threshold or higher)
  datatest=data[exp(data$BTTLDtotLN)>i,]
  # Use zelig to estimate hazard 
  z.out=zelig(Surv(Duration,ongoingCONF)~DRUGP+mountain+resources+gdplLN+
                Weakrebel+Equalstrenght+internatz+coldwar+popLN+ethfrac+DEM+
                terr+confbord,model="coxph",data=datatest,robust=T,cluster="ccode")
  # Set treatment and control
  x.high=setx(z.out,DRUGP=1)
  x.low=setx(z.out,DRUGP=0)
  # Simulate hazards for treated and control data
  set.seed(4321)
  s.out=sim(z.out,x=x.low,x1=x.high)
  # Extract percent change in hazards as well as confidence intervals
  fd.est=c(fd.est,mean(s.out$qi$hr))
  low.est=c(low.est,quantile(s.out$qi$hr,.025))
  high.est=c(high.est,quantile(s.out$qi$hr,.975))
  low90.est=c(low90.est,quantile(s.out$qi$hr,.05))
  high90.est=c(high90.est,quantile(s.out$qi$hr,.95))
  nd.est=c(nd.est,dim(datatest)[1])
  print(i)
}

# Create the figure by plotting the mean percent changes and confidence intervals
par(mar=c(5,5,2,2)) 
par(mgp=c(2.2,1,0))
plot(btldt,fd.est,ylab="Simulated percent Change in the Hazard of Conflict Termination\n(No Drug Cultivation in Conflict Zone to Drug Cultivation in Conflict Zone)",ylim=c(round(min(low.est),1),1.6),
     ,pch=20,xlab="Battle Death Threshold",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)
lines(btldt,low.est,lty=1,col="gray39",lwd=2)
lines(btldt,high.est,lty=1,col="gray39",lwd=2)
lines(btldt,low90.est,lty=1,col="gray",lwd=2)
lines(btldt,high90.est,lty=1,col="gray",lwd=2)
axis(1, at=seq(0,1000,50), labels=seq(0,1000,50),cex.axis=.7)
axis(2, at=seq(.1,1.6,.1), 
     labels=seq(-90,60,10),cex.axis=.65)
abline(h=1,lty=1)
legend(x=620,y=1.55,legend=c("95% Confidence Interval    ","90% Confidence Interval"),
       lty=c(1,1),lwd=c(1.8,1.8), col=c("gray39","gray"),cex=.66)

#########################################
##### FIGURE 6 (Conflict Intensity) #####
#########################################

### Part A: Excluding duration control 

# Select variables for the analysis and pare the dataset down
keep.vars <- c("BTTLDintLN","ALLGEMP","DRUGP","popLN","hydrocarbonP",
               "OUTZhydrocarbonP","mountain","durationLN","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","BTTLDtotLN","ethPOL","DEM", "ccode")  
all.data <- dataset[keep.vars]
data <- na.omit(all.data)

# Create empty objects to record percent changes and confidence intervals
fdnc.estint=c()
lownc.estint=c()
highnc.estint=c()
low90nc.estint=c()
high90nc.estint=c()
n.estintnc=c()
t.estintnc=c()

for (i in btldt){
  # Create dataset for the new battle death threshold (i.e. keep data with only the threshold or higher)
  datatest=data[exp(data$BTTLDtotLN)>i,]
  # Create X covariates for the new dataset
  X=datatest[,c("ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord",
                "Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")]
  # Get percent change and confidence intervals going from control to treated by using the qsim function created above.   
  set.seed(4321)
  tester=qsim(X,datatest,treat)
  fdnc.estint=c(fdnc.estint,tester[1])
  lownc.estint=c(lownc.estint,tester[2])
  highnc.estint=c(highnc.estint,tester[3])
  low90nc.estint=c(low90nc.estint,tester[4])
  high90nc.estint=c(high90nc.estint,tester[5])
  n.estintnc=c(n.estintnc,dim(datatest)[1])
  t.estintnc=c(t.estintnc,sum(datatest$DRUGP))
  print(i)
}

### Part B: Including duration control

fd.est=c()
low.est=c()
high.est=c()
low90.est=c()
high90.est=c()
n.est=c()
t.est=c()

for (i in btldt){
  # Create dataset for the new battle death threshold (i.e. keep data with only the threshold or higher) 
  datatest=data[exp(data$BTTLDtotLN)>i,]
  # Select independent variables
  X=datatest[,c("ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM","durationLN")]
  # Get percent change and confidence intervals going from control to treated by using the qsim function created above. 
  set.seed(4321)
  tester=qsim(X,datatest,treat)
  fd.est=c(fd.est,tester[1])
  low.est=c(low.est,tester[2])
  high.est=c(high.est,tester[3])
  low90.est=c(low90.est,tester[4])
  high90.est=c(high90.est,tester[5])
  n.est=c(n.est,dim(datatest)[1])
  t.est=c(t.est,sum(datatest$treat))
  print(i)
}

# Create Figure 6a by plotting the mean percent changes and confidence intervals
par(mar=c(5,5,2,2)) 
par(mgp=c(2.2,1,0))
plot(btldt,fdnc.estint,ylim=c(round(min(lownc.estint),1),round(max(high.est),1)),ylab="Percent Change in in Battle Deaths per Year\n(No Drug Cultivation in Conflict Zone to Drug Cultivation in Conflict Zone)"
     ,pch=20,xlab="Battle Death Threshold",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)
lines(btldt,lownc.estint,lty=1,col="gray39",lwd=2)
lines(btldt,highnc.estint,lty=1,col="gray39",lwd=2)
lines(btldt,low90nc.estint,lty=1,col="gray",lwd=2)
lines(btldt,high90nc.estint,lty=1,col="gray",lwd=2)
axis(1, at=seq(0,1000,50), labels=seq(0,1000,50),cex.axis=.7)
axis(2, at=seq(-2600,1400,200), labels=seq(-2600,1400,200),cex.axis=.7)
axis(2, at=seq(round(min(lownc.estint),1),round(max(high.est),1),.1),labels=seq(round(min(lownc.estint),1)*100,round(max(high.est),1)*100,.1*100), cex.axis=.7)
abline(h=0,lty=1)
legend(x=10,y=.63,legend=c("95% Confidence Interval    ","90% Confidence Interval"),
       lty=c(1,1),lwd=c(1.8,1.8), col=c("gray39","gray"),cex=.5)

# Create Figure 6b by plotting the mean percent changes and confidence intervals
par(mar=c(5,5,2,2)) 
par(mgp=c(2.2,1,0))
plot(btldt,fd.est,ylim=c(round(min(lownc.estint),1),round(max(high.est),1)),ylab="Percent Change in in Battle Deaths per Year\n(No Drug Cultivation in Conflict Zone to Drug Cultivation in Conflict Zone)",pch=20,xlab="Battle Death Threshold",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)
lines(btldt,low.est,lty=1,col="gray39",lwd=2)
lines(btldt,high.est,lty=1,col="gray39",lwd=2)
lines(btldt,low90.est,lty=1,col="gray",lwd=2)
lines(btldt,high90.est,lty=1,col="gray",lwd=2)
axis(1, at=seq(0,1000,50), labels=seq(0,1000,50),cex.axis=.7)
axis(2, at=seq(round(min(lownc.estint),1),round(max(high.est),1),.1),labels=seq(round(min(lownc.estint),1)*100,round(max(high.est),1)*100,.1*100), cex.axis=.7)
abline(h=0,lty=1)
legend(x=10,y=.63,legend=c("95% Confidence Interval    ","90% Confidence Interval"),lty=c(1,1),lwd=c(1.8,1.8), col=c("gray39","gray"),cex=.5)

###################################################################
##### APPENDIX E (Random Drop Tests for Duration & Intensity) #####
###################################################################

## APPENDIX E.1: Conflict Duration

# Select variables for the analysis and pare the dataset down
keep.vars <- c("BTTLDtotLN","BTTLDtot","resources","ongoingCONF","Duration","DRUGP","mountain","confbord","Weakrebel","Equalstrenght","coldwar",
               "internatz","popLN","ethfrac","DEM","ccode","terr","gdplLN")  
all.data <- dataset[keep.vars]
data <- na.omit(all.data)

# Create censored observations
data$ongoingCONF=1-data$ongoingCONF

# Simulate the effect of simply randomly deleting as many observations as are deleted at each threshold above 
# by sampling 1000 datasets for each threshold

# First, find number of values deleted at each threshold
index=c()
for (i in 2:length(nd.est)-1){
  index[i]=(dim(data)[1]-nd.est[i])
}
deld.est=unique(index)

# Create empty objects for the percent changes and confidence intervals
fddd.rans=c()
lowdd.rans=c()
highdd.rans=c()
lowdd.rans90=c()
highdd.rans90=c()
nd.rans=c()
td.rans=c()

hold=cbind(btldt,nd.est)
hold=rbind(1,hold)
x=matrix(,dim(hold)[1],3)
for (i in 2:dim(hold)[1]){
  x[i,1:3]=ifelse(hold[i,2]==hold[i-1,2],NA,i)
}

temporary=na.omit(x[,1])
uniquenh=(hold[temporary,])

# Start from the 2nd observation as the first is zero
for (i in 2:length(deld.est)){
  fd.ransin=c()
  low.ransin=c()
  high.ransin=c()
  lowransin90.est=c()
  highransin90.est=c()
  n.ransin=c()
  t.ransin=c()
  # Now create 1000 samples of deleted observations
  for (j in 1:1000){
    # Randomly select the number of data points deleted for each threshold
    set.seed(j+(i*1000))
    rand=sample(dim(data)[1],deld.est[i], replace=F) 
    # Remove randomly selected observations from the dataset
    datatest=data[-rand,]
    # As above, use zelig to simulate percent changes
    z.out=zelig(Surv(Duration,ongoingCONF)~DRUGP+mountain+resources+gdplLN+Weakrebel+Equalstrenght+internatz+
                  coldwar+popLN+ethfrac+DEM+terr+confbord,model="coxph",data=datatest,robust=T,cluster="ccode")   
    x.high=setx(z.out,DRUGP=1)
    x.low=setx(z.out,DRUGP=0)
    s.out=sim(z.out,x=x.low,x1=x.high)
    fd.ransin=c(fd.ransin,mean(s.out$qi$hr))
    low.ransin=c(low.ransin,quantile(s.out$qi$hr,.025))
    high.ransin=c(high.ransin,quantile(s.out$qi$hr,.975))
    lowransin90.est=c(lowransin90.est,quantile(s.out$qi$hr,.05))
    highransin90.est=c(highransin90.est,quantile(s.out$qi$hr,.95))
  }
  fddd.rans=c(fddd.rans,mean(fd.ransin))
  lowdd.rans=c(lowdd.rans,mean(low.ransin))
  highdd.rans=c(highdd.rans,mean(high.ransin))
  lowdd.rans90=c(lowdd.rans90,mean(lowransin90.est))
  highdd.rans90=c(highdd.rans90,mean(highransin90.est))
  print(i)
}

# Create the plot by graphing the mean percentage change across each of the 1,000 datasets at each drop threshold
par(mar=c(5,5,4,2),mgp=c(2.2,1,0))
plot(1:28,fddd.rans,ylab="Percent Change in the Hazard of Conflict Termination\n(No Drugs in Conflict Zone to Drugs in Conflict Zone)",ylim=c(.2,2.35),
     pch=20,xlab="Number of Observations Randomly Deleted",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)
lines(1:28,lowdd.rans,lty=1,col="gray39",lwd=2)
lines(1:28,highdd.rans,lty=1,col="gray39",lwd=2)
lines(1:28,lowdd.rans90,lty=1,col="gray",lwd=2)
lines(1:28,highdd.rans90,lty=1,col="gray",lwd=2)
axis(1, at=seq(1,28,1), labels=F,cex.axis=.7,tick=T)
axis(1, at=seq(1,28,2), labels=deld.est[ seq(2,28,2)],cex.axis=.7,tick=F)
axis(2, at=seq(.4,2.8,.1), labels=seq(-60,180,10),cex.axis=.65)
axis(3, at=1:28,label=F,cex.axis=.7)
axis(3, at=seq(1,28,2),label=uniquenh[seq(2,28,2),1],cex.axis=.7,tick=F)
mtext("Corresponding Battle Death Threshold", side=3, cex=.8,line=2)
abline(h=0,lty=1)
abline(h=1,lty=1)
legend(x=17,y=2.78,legend=c("95% Confidence Interval    ","90% Confidence Interval"),
       lty=c(1,1),lwd=c(1.8,1.8), col=c("gray39","gray"),cex=.66)

## APPENDIX E.2: Conflict Intensity

# Select variables for the analysis and pare the dataset down
keep.vars <- c("BTTLDintLN","ALLGEMP","DRUGP","popLN","hydrocarbonP",
               "OUTZhydrocarbonP","mountain","durationLN","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","BTTLDtotLN","ethPOL","DEM", "ccode")  
all.data <- dataset[keep.vars]
data <- na.omit(all.data)

## Figure E.2a: excluding duration control

# Simulate the effect of simply randomly deleting as many observations as are deleted at each threshold above 
# by sampling 1000 datasets for each threshold

# First, find number of values deleted at each threshold
index=c()
for (i in 2:length(n.est)-1){
  index[i]=(dim(data)[1]-n.est[i])
}

# Get unique set of number of deleted observations
deld.est=unique(index)

# Create empty objects for the percent changes and confidence intervals
fddnd.rans=c()
lowdnd.rans=c()
highdnd.rans=c()
lowdnd.rans90=c()
highdnd.rans90=c()
nnd.rans=c()
tnd.rans=c()

# Start from the 2nd observation as the first is zero
for (i in 2:length(deld.est)){
  # Create empty objects for each of the 1,000 datasets to be produced at each drop value  
  fd.ransin=c()
  low.ransin=c()
  high.ransin=c()
  low.ransin90=c()
  high.ransin90=c()
  n.ransin=c()
  t.ransin=c()
  # Now create 1000 samples of deleted observations
  for (j in 1:1000){
    # Randomly select the number of data points to be deleted
    set.seed(j+(i*1000))
    rand=sample(dim(data)[1],deld.est[i], replace=F) 
    # Remove randomly selected observations from the dataset
    datatest=data[-rand,]
    X=datatest[,c("ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")]
    tester=qsim(X,datatest,treat)
    fd.ransin=c(fd.ransin,tester[1])
    low.ransin=c(low.ransin,tester[2])
    high.ransin=c(high.ransin,tester[3])
    low.ransin90=c(low.ransin90,tester[4])
    high.ransin90=c(high.ransin90,tester[5])
    n.ransin=c(n.ransin,dim(datatest)[1])
    t.ransin=c(t.ransin,sum(datatest$DRUGP))
  }
  fddnd.rans=c(fddnd.rans,mean(fd.ransin))
  lowdnd.rans=c(lowdnd.rans,mean(low.ransin))
  highdnd.rans=c(highdnd.rans,mean(high.ransin))
  lowdnd.rans90=c(lowdnd.rans90,mean(low.ransin90))
  highdnd.rans90=c(highdnd.rans90,mean(high.ransin90))
  nnd.rans=c(nnd.rans,mean(n.ransin))
  tnd.rans=c(tnd.rans,round(mean(t.ransin),0))
  print(i)
}

summary=cbind(btldt,nd.est,t.est)
summary=rbind(1,summary)
x=matrix(,dim(summary)[1],3)
for (i in 2:dim(summary)[1]){
  x[i,1:3]=ifelse(summary[i,2]==summary[i-1,2],NA,i)
}
temporary=na.omit(x[,1])
uniquend=(summary[temporary,])

# Create the plot by graphing the mean percentage change across each of the 1,000 datasets at each drop threshold
par(mar=c(5,5,4,2),mgp=c(2.5,1,0))
par(mgp=c(2.2,1,0))
plot(1:28,fddnd.rans,ylim=c(-.92,.6),ylab="Percent Change in in Battle Deaths per Year\n(No Drug Cultivation in Conflict Zone to Drug Cultivation in Conflict Zone)",
     pch=20,xlab="Battle Death Threshold",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)
lines(1:28,lowdnd.rans,lty=1,col="gray39",lwd=2)
lines(1:28,highdnd.rans,lty=1,col="gray39",lwd=2)
lines(1:28,lowdnd.rans90,lty=1,col="gray",lwd=2)
lines(1:28,highdnd.rans90,lty=1,col="gray",lwd=2)
axis(1, at=seq(1,28,1), labels=F,cex.axis=.7,tick=T)
axis(1, at=seq(1,28,2), labels=deld.est[ seq(2,28,2)],cex.axis=.7,tick=F)
axis(2, at=seq(-.9,.6,.1),labels=seq(-90,60,10), cex.axis=.7)
axis(3, at=1:28,label=F,cex.axis=.7)
axis(3, at=seq(1,28,2),label=uniquend[seq(2,28,2),1],cex.axis=.7,tick=F)
mtext("Corresponding Battle Death Threshold", side=3, cex=.8,line=2)
abline(h=0,lty=1)
abline(h=0,lty=1)
legend(x=.5,y=.6,legend=c("95% Confidence Interval    ","90% Confidence Interval"),
       lty=c(1,1),lwd=c(1.8,1.8), col=c("gray39","gray"),cex=.66)

## Figure E.2b: including duration control

# Simulate the effect of simply randomly deleting as many observations as are deleted at each threshold above 
# by sampling 1000 datasets for each threshold

# First, find number of values deleted at each threshold
index=c()
for (i in 2:length(n.est)-1){
  index[i]=(dim(data)[1]-n.est[i])
}

# Get unique set of number of deleted observations
deld.est=unique(index)

# Create empty objects for the percent changes and confidence intervals
fdd.rans=c()
lowd.rans=c()
highd.rans=c()
lowd.rans90=c()
highd.rans90=c()
n.rans=c()
t.rans=c()

# Start from the 2nd observation as the first is zero
for (i in 2:length(deld.est)){
  # Create empty objects for each of the 1,000 datasets to be produced at each drop value
  fd.ransin=c()
  low.ransin=c()
  high.ransin=c()
  low.ransin90=c()
  high.ransin90=c()
  n.ransin=c()
  t.ransin=c()
  # Now create 1000 samples of deleted observations
  for (j in 1:1000){
    # Randomly select the number of data points to be deleted
    set.seed(j+(i*1000))
    rand=sample(dim(data)[1],deld.est[i], replace=F) 
    # Remove randomly selected observations from the dataset
    datatest=data[-rand,]
    X=datatest[,c("ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM","durationLN")]
    tester=qsim(X,datatest,treat)
    fd.ransin=c(fd.ransin,tester[1])
    low.ransin=c(low.ransin,tester[2])
    high.ransin=c(high.ransin,tester[3])
    low.ransin90=c(low.ransin90,tester[4])
    high.ransin90=c(high.ransin90,tester[5])
    n.ransin=c(n.ransin,dim(datatest)[1])
    t.ransin=c(t.ransin,sum(datatest$DRUGP))
  }
  fdd.rans=c(fdd.rans,mean(fd.ransin))
  lowd.rans=c(lowd.rans,mean(low.ransin))
  highd.rans=c(highd.rans,mean(high.ransin))
  lowd.rans90=c(lowd.rans90,mean(low.ransin90))
  highd.rans90=c(highd.rans90,mean(high.ransin90))
  n.rans=c(n.rans,mean(n.ransin))
  t.rans=c(t.rans,round(mean(t.ransin),0))
  print(i)
}

summary=cbind(btldt,nd.est,t.est)
summary=rbind(1,summary)
x=matrix(,dim(summary)[1],3)
for (i in 2:dim(summary)[1]){
  x[i,1:3]=ifelse(summary[i,2]==summary[i-1,2],NA,i)
}

# Get unique number of deleted observations
temporary=na.omit(x[,1])
uniquen=(summary[temporary,])

# Create the plot by graphing the mean percentage change across each of the 1,000 datasets at each drop threshold
par(mar=c(5,5,4,2),mgp=c(2.2,1,0))
plot(1:28,fdd.rans,ylim=c(-.92,.62),ylab="Percent Change in in Battle Deaths per Year\n(No Drug Cultivation in Conflict Zone to Drug Cultivation in Conflict Zone)"
     ,pch=20,xlab="Battle Death Threshold",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)
lines(1:28,lowd.rans,lty=1,col="gray39",lwd=2)
lines(1:28,highd.rans,lty=1,col="gray39",lwd=2)
lines(1:28,lowd.rans90,lty=1,col="gray",lwd=2)
lines(1:28,highd.rans90,lty=1,col="gray",lwd=2)
axis(1, at=seq(1,28,1), labels=F,cex.axis=.7,tick=T)
axis(1, at=seq(1,28,2), labels=deld.est[ seq(2,28,2)],cex.axis=.7,tick=F)
axis(2, at=seq(-.9,.6,.1),labels=seq(-90,60,10), cex.axis=.7)
axis(3, at=1:28,label=F,cex.axis=.7)
axis(3, at=seq(1,28,2),label=uniquen[seq(2,28,2),1],cex.axis=.7,tick=F)
mtext("Corresponding Battle Death Threshold", side=3, cex=.8,line=2)
abline(h=0,lty=1)
abline(h=0,lty=1)
legend(x=.5,y=.6,legend=c("95% Confidence Interval    ","90% Confidence Interval"),
       lty=c(1,1),lwd=c(1.8,1.8), col=c("gray39","gray"),cex=.66)

##########################################
##### APPENDIX F (Matching Analysis) #####
##########################################

# Replicate function from above, but add in option for whether or not regression should include weighting or
# complete regressions with clustered standard errors and simulate percent change in dependent variable. 
# Function is composed of three elements: the independent variables, the dataset, and indicator for the treatment variable.  
# Note that the DV should be the first column in the data given to the function.

qsim=function(X,data,treat,fact){
  # Regress data and given covariates
  if(fact==0){M2b <- lm(data[,1]~as.matrix(X), data=data) }
  if(fact==1){M2b <- lm(data[,1]~as.matrix(X),weight=weights, data=data) }
  # Generate a robust covariance matrix, clustered by country
  rvarcov=clv(data,M2b,data$ccode)
  # Create design matrix
  Xdes=as.matrix(cbind(1,X))
  # Create covariates values to simulate percent change
  covs=as.matrix(apply(Xdes,2,mean))
  # Create "treatment"
  covs[treat,]<-1
  # Create "control"
  covs2=covs
  covs2[treat,]<-0
  # Create empty object for the treated dataset
  simdata=c()
  # Create empty object for the control dataset
  simdata2=c()
  for(i in 1:1000){
    # Simulate betas from the original model
    simbetas=rmvnorm(1,mean=M2b$coefficients,sigma=rvarcov)
    simsigma=sqrt(M2b$df.residual* summary(M2b)$sigma^2 / rchisq(1, M2b$df.residual))
    # Generate the simulated percent change values
    simdata[i]=(mean(rnorm(n=1000,(simbetas%*%covs),(simsigma))))
    simdata2[i]=(mean(rnorm(n=1000,(simbetas%*%covs2),(simsigma))))
  }
  # Create vector of percent change 
  TEs=(exp(simdata)-exp(simdata2))/exp(simdata2)
  # Take average of the 1,000 percent change values
  ate=(mean(TEs))
  # Confidence intervals
  low=(quantile(TEs,.025))              
  high=(quantile(TEs,.975))  
  low90=(quantile(TEs,.05))              
  high90=(quantile(TEs,.95)) 
  # Print the percent change and confidence intervals
  print(cbind(ate,low,high,low90,high90))
}

## TOTAL BATTLE DEATHS: Figure F.3a

# Select variables for the analysis and pare the dataset down
keep.vars <- c("BTTLDtotLN","BTTLDtot","ccode","DRUGP","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")
all.data <- dataset[keep.vars]
data <- na.omit(all.data)

# Create empty objects to record percent changes and confidence intervals
fd.est=c()
low.est=c()
high.est=c()
low90.est=c()
high90.est=c()

# Run the analysis with no matching and also with various matching algorithms
method=c("none","nearest","optimal","full")
for (i in 1:length(method)){
  set.seed(4321)
  if (i>1) {matched <- matchit(DRUGP~ALLGEMP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+
                                                 Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+
                                                 DEM, data=data, method=method[i])}
  if (i==1) {datatest=data}
   if (i>1) {datatest=match.data(matched)}
  X=datatest[,c("ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")]
  # Get percent change and confidence intervals going from control to treated by using the qsim function created above
  set.seed(4321)
  if (i<3) { tester=qsim(X,datatest,treat,0) }
  if (i>2) { tester=qsim(X,datatest,treat,1) }
  fd.est=c(fd.est,tester[1])
  low.est=c(low.est,tester[2])
  high.est=c(high.est,tester[3])
  low90.est=c(low90.est,tester[4])
  high90.est=c(high90.est,tester[5])
  print(i)  
}

# Plot results
par(mar=c(3,5,2,2)) 
par(mgp=c(2.2,1,0))
plot(seq(1:5),seq(1:5),xlim=c(.5,5.5),bty="n",ylim=c(min(low.est),.3),ylab="Percent Change in Battle Deaths \n(No Drug Cultivation in Conflict Zone to Drug Cultivation in Conflict Zone)",
     pch=20,cex.lab=.7,cex.axis=.7,yaxt = "n",xaxt = "n",type='n',xlab="")
arrows(c(1,3,4,5),low.est[c(1,2,3,4)],c(1,3,4,5),high.est[c(1,2,3,4)],lty=1,col="gray39",lwd=2.4, angle=90, code=3,length = 0.05)
arrows(c(1,3,4,5),low90.est[c(1,2,3,4)],c(1,3,4,5),high90.est[c(1,2,3,4)],lty=1,col="gray",lwd=2.4, angle=90, code=3,length = 0.05)
points(c(1,3,4,5),fd.est[c(1,2,3,4)],pch=20,cex=1.2)
segments(.5,0,5.5,0,lty=1)
legend("topright",legend=c("95% Confidence Interval    ","90% Confidence Interval"),lty=c(1,1),lwd=c(2,2), col=c("gray39","gray"),cex=.66)
axis(2, at=seq(round(min(low.est),1),.3,.1),cex.axis=.7,labels=seq(100*round(min(low.est),1),100*.3,100*.1))
mtext("Un-Matched", at=par("usr")[1]+.7,side=1, line=-.15, cex.lab=.72,las=1,cex=.72,padj = 1)
mtext(substitute(paste("(",italic('N'), " = 244)" )), at=par("usr")[1]+.7,side=1, line=.55, cex.lab=.72,las=1,cex=.72,padj = 1)
axis(side = 1,at=c(.7,1.28),labels = F,tcl=0,line=.45,lwd=1.2)
mtext("Nearest", at=par("usr")[1]+2.7,side=1, line=-.15, cex.lab=.72,las=1,cex=.72,padj = 1)
mtext(substitute(paste("(",italic('N'), " = 68)" )), at=par("usr")[1]+2.7,side=1, line=.55, cex.lab=.72,las=1,cex=.72,padj = 1)
axis(side = 1,lwd=1.2,at=c(2.7,3.28),labels = F,tcl=0,line=.45)
mtext("Optimal", at=par("usr")[1]+3.7,side=1, line=-.15, cex.lab=.72,las=1,cex=.72,padj = 1)
mtext(substitute(paste("(",italic('N'), " = 68)" )), at=par("usr")[1]+3.7,side=1, line=.55, cex.lab=.72,las=1,cex=.72,padj = 1)
axis(side = 1,lwd=1.2,at=c(3.7,4.28),labels = F,tcl=0,line=.45)
mtext("Full", at=par("usr")[1]+4.7,side=1, line=-.15, cex.lab=.72,las=1,cex=.72,padj = 1)
mtext(substitute(paste("(",italic('N'), " = 244)" )), at=par("usr")[1]+4.7,side=1, line=.55, cex.lab=.72,las=1,cex=.72,padj = 1)
axis(side = 1,at=c(4.7,5.28),labels = F,tcl=0,line=.45,lwd=1.2)

## DURATION: Figure F.3b

# Select variables for the analysis and pare the dataset down
keep.vars <- c("BTTLDtotLN","BTTLDtot","resources","ongoingCONF","Duration","DRUGP","mountain","confbord","Weakrebel","Equalstrenght","coldwar",
               "internatz","popLN","ethfrac","DEM","ccode","terr","gdplLN")  
data <- na.omit(dataset[keep.vars])

# Create censored observations
data$ongoingCONF=1-data$ongoingCONF

# Create empty objects to record percent changes and confidence intervals
fd.est=c()
low.est=c()
high.est=c()
low90.est=c()
high90.est=c()
n.est=c()
t.est=c()

# Run the analysis with no matching and also with various matching algorithms
method=c("none","nearest","optimal","full")
for (i in 1:length(method)){
  set.seed(4321)
  if (i>1) {matched <- matchit(DRUGP~mountain+resources+gdplLN+Weakrebel+Equalstrenght+internatz+coldwar+popLN+
                                                 ethfrac+DEM+terr+confbord, data=data, method=method[i])}
  if (i==1) {datatest=data}
  if (i>1) {datatest=match.data(matched)}
  if (i>1) {matweight=as.vector(datatest$weights)}
  if (i==1) {z.out=zelig(Surv(Duration,ongoingCONF)~DRUGP+mountain+resources+gdplLN+Weakrebel+Equalstrenght+internatz+coldwar+popLN+ethfrac+DEM+terr+confbord,
                        model="coxph",data=datatest,robust=T,cluster="ccode")}
  if (i>3) {z.out=zelig(Surv(Duration,ongoingCONF)~DRUGP+mountain+resources+gdplLN+Weakrebel+Equalstrenght+internatz+coldwar+popLN+ethfrac+DEM+terr+confbord,
                        model="coxph",data=datatest,weights=matweight,robust=T,cluster="ccode")}
  if (i==2|i==3) {z.out=zelig(Surv(Duration,ongoingCONF)~DRUGP+mountain+resources+gdplLN+Weakrebel+Equalstrenght+internatz+coldwar+popLN+ethfrac+DEM+terr+confbord,
                        model="coxph",data=datatest,robust=T,cluster="ccode")}
  x.high=setx(z.out,DRUGP=1)
  x.low=setx(z.out,DRUGP=0)
  set.seed(4321)
  s.out=sim(z.out,x=x.low,x1=x.high)
  fd.est=c(fd.est,mean(s.out$qi$hr))
  low.est=c(low.est,quantile(s.out$qi$hr,.025))
  high.est=c(high.est,quantile(s.out$qi$hr,.975))
  low90.est=c(low90.est,quantile(s.out$qi$hr,.05))
  high90.est=c(high90.est,quantile(s.out$qi$hr,.95))
  n.est=c(n.est,dim(datatest)[1])
}

# Plot results
par(mar=c(3,5,2,2)) 
par(mgp=c(2.2,1,0))
plot(seq(1:5),seq(1:5),xlim=c(.5,5.5),bty="n",ylim=c(min(low.est),1.55),ylab="Percent Change in the Hazard of Conflict Termination\n(No Drug Cultivation in Conflict Zone to Drug Cultivation in Conflict Zone)",pch=20,cex.lab=.7,cex.axis=.7,yaxt = "n",xaxt = "n",type='n',xlab="")
segments(.5,1,5.5,1,lty=1)
arrows(c(1,3,4,5),low.est[c(1,2,3,4)],c(1,3,4,5),high.est[c(1,2,3,4)],lty=1,col="gray39",lwd=2.4, angle=90, code=3,length = 0.05)
arrows(c(1,3,4,5),low90.est[c(1,2,3,4)],c(1,3,4,5),high90.est[c(1,2,3,4)],lty=1,col="gray",lwd=2.4, angle=90, code=3,length = 0.05)
points(c(1,3,4,5),fd.est[c(1,2,3,4)],pch=20,cex=1.2)
legend("topright",legend=c("95% Confidence Interval","90% Confidence Interval"),lty=c(1,1),lwd=c(2,2), col=c("gray39","gray"),cex=.5)
axis(2, at=seq(.3,1.5,.1), labels=seq(-70,50,10),cex.axis=.65)
mtext("Un-Matched", at=par("usr")[1]+.7,side=1, line=-.15, cex.lab=.72,las=1,cex=.72,padj = 1)
mtext(substitute(paste("(",italic('N'), " = 238)" )), at=par("usr")[1]+.7,side=1, line=.55, cex.lab=.72,las=1,cex=.72,padj = 1)
axis(side = 1,at=c(.7,1.28),labels = F,tcl=0,line=.45,lwd=1.2)
mtext("Nearest", at=par("usr")[1]+2.7,side=1, line=-.15, cex.lab=.72,las=1,cex=.72,padj = 1)
mtext(substitute(paste("(",italic('N'), " = 68)" )), at=par("usr")[1]+2.7,side=1, line=.55, cex.lab=.72,las=1,cex=.72,padj = 1)
axis(side = 1,lwd=1.2,at=c(2.7,3.28),labels = F,tcl=0,line=.45)
mtext("Optimal", at=par("usr")[1]+3.7,side=1, line=-.15, cex.lab=.72,las=1,cex=.72,padj = 1)
mtext(substitute(paste("(",italic('N'), " = 68)" )), at=par("usr")[1]+3.7,side=1, line=.55, cex.lab=.72,las=1,cex=.72,padj = 1)
axis(side = 1,lwd=1.2,at=c(3.7,4.28),labels = F,tcl=0,line=.45)
mtext("Full", at=par("usr")[1]+4.7,side=1, line=-.15, cex.lab=.72,las=1,cex=.72,padj = 1)
mtext(substitute(paste("(",italic('N'), " = 238)" )), at=par("usr")[1]+4.7,side=1, line=.55, cex.lab=.72,las=1,cex=.72,padj = 1)
axis(side = 1,at=c(4.7,5.28),labels = F,tcl=0,line=.45,lwd=1.2)